function [HHstates,wage_moments,transitions]=data_simulation_A12_end(Par1,ff_2,wf_2,fi_2,wi_2,Wf,Wi,Wn)


df_2=Par1(1);
di_2=Par1(2);
lnf_2=Par1(3);
lni_2=Par1(4);
lfi_2=Par1(5);
lif_2=Par1(6);
alphaf=Par1(7);
betaf=Par1(8);
alphai=Par1(9);
betai=Par1(10);
% Data simulation (MSM)

N_households=10000; % # households
maxtime=400;   % # quarters
rng('default'); % set seed

% Pre-allocation
value=zeros(N_households,maxtime);
state=zeros(N_households,maxtime);
wage_draw=ones(N_households,maxtime);

for iter=1:N_households
    
    time=1; % households are born in the situation 1: "N"
    scale = lni_2  + lnf_2  + (1 - lni_2)*(1-lnf_2);
    randnumber=scale*rand;
    
    if randnumber < lnf_2
        urn = random('Beta', alphaf, betaf);
        transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
        wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
        value(iter,time)=Wf(wage_aux)*(Wf(wage_aux)>Wn)+Wn*(Wf(wage_aux)<=Wn);
        state(iter,time)=(Wf(wage_aux)<=Wn)+2*(Wf(wage_aux)>Wn);
        wage_draw(iter,time) = wage_aux*(Wf(wage_aux)>Wn) + (Wf(wage_aux)<=Wn);
        
    elseif randnumber >=lnf_2 && randnumber<(lni_2+lnf_2)
        urn = random ('Beta', alphai, betai);
        transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
        wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
        value(iter,time)=Wi(wage_aux)*(Wi(wage_aux)>Wn)+Wn*(Wi(wage_aux)<=Wn);
        state(iter,time)=3*(Wi(wage_aux)>Wn)+ (Wi(wage_aux)<=Wn);
        wage_draw(iter,time) = wage_aux*(Wi(wage_aux)>Wn) + (Wi(wage_aux)<=Wn);
         
    else
        wage_draw(iter,time)=1;
        value(iter,time)=Wn;
        state(iter,time)=1;
        
    end
    
    
    time=2;
    
    while time<=maxtime
        
        
        
        
        %% F in t-1
        if (state(iter,time-1) == 2)
            scale = df_2 + lfi_2 + (1 - df_2)*(1 - lfi_2);
            randnumber=scale*rand;
            if  (randnumber < df_2)
                wage_draw(iter,time)=1;
                value(iter,time)=Wn;
                state(iter,time)=1;
            end
            if randnumber >= df_2 && randnumber<(lfi_2+df_2)

                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                value(iter,time)=Wi(wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=3*(value(iter,time)> value(iter,time-1))+2*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw(iter,time) = wage_draw(iter,time - 1);
                else
                    wage_draw(iter,time) = wage_aux;
                end
            end
            if randnumber>= (lfi_2+df_2)
  		wage_draw(iter,time)=wage_draw(iter,time-1);
                value(iter,time)=value(iter,time-1);
                state(iter,time)=state(iter,time-1);
            end
        end
        %% I em t-1
        if state(iter,time-1)==3
            scale = di_2 + lif_2 + (1-di_2)*(1-lif_2);
            randnumber=scale*rand;
            
            if randnumber < di_2
                wage_draw(iter,time)=1;
                value(iter,time)=Wn;
                state(iter,time)=1;
            end
   
            if randnumber >= (di_2) && randnumber<(lif_2+di_2)
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
                wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
                value(iter,time)=Wf(wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=2*(value(iter,time)> value(iter,time-1))+3*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw(iter,time) = wage_draw(iter,time - 1);
                else
                    wage_draw(iter,time) = wage_aux;
                end
            end
            if randnumber>=(lif_2+di_2)
                wage_draw(iter,time)=wage_draw(iter,time-1);
                value(iter,time)=value(iter,time-1);
                state(iter,time)=state(iter,time-1);
            end
        end
       

        % N in t-1
        if state(iter,time-1) ==1
            scale = lnf_2 + lni_2 +(1 - lni_2)*(1 - lnf_2);
            randnumber=scale*rand;
            
            if randnumber < lnf_2
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
                wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
                value(iter,time)=Wf(wage_aux)*(Wf(wage_aux)>Wn)+Wn*(Wf(wage_aux)<=Wn);
                state(iter,time)=2*(Wf(wage_aux)>Wn)+(Wf(wage_aux)<=Wn);
                 if state (iter,time) == state(iter,time-1)
                    wage_draw(iter,time) = wage_draw(iter,time - 1);
                else
                    wage_draw(iter,time) = wage_aux;
                end
            end
            if randnumber >=lnf_2 && randnumber<(lni_2+lnf_2)
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                value(iter,time)=Wi(wage_aux)*(Wi(wage_aux)>Wn)+Wn*(Wi(wage_aux)<=Wn);
                state(iter,time)=3*(Wi(wage_aux)>Wn)+ (Wi(wage_aux)<=Wn);
                 if state (iter,time) == state(iter,time-1)
                    wage_draw(iter,time) = wage_draw(iter,time - 1);
                else
                    wage_draw(iter,time) = wage_aux;
                end
            end
            
            if randnumber>=(lni_2+lnf_2)
                value(iter,time)=Wn;
                state(iter,time)=1;
            end
        end
        
        
        time=time+1;
        
    end
    
end

% Moments calculation

% use only last 2 quarters of simulation for each household
NF=sum(state(:,end-1)==2)/N_households;
NI=sum(state(:,end-1)==3)/N_households;
NN=sum(state(:,end-1)==1)/N_households;


HHstates=[NF;NI;NN]; % joint-job state

wage_formal_2=wf_2(wage_draw(:,end-1));
wage_informal_2=wi_2(wage_draw(:,end-1));


indicewf_2=find(state(:,end-1)==2);
wage_formal_2=sort(wage_formal_2(indicewf_2));
indicewi_2=find(state(:,end-1)==3);
wage_informal_2=sort(wage_informal_2(indicewi_2));

sdwi_2 = std(wage_informal_2);

sdwf_2 = std(wage_formal_2);


meanwf_2=mean(wage_formal_2);

p10wf_2=prctile(wage_formal_2,10);
p25wf_2=prctile(wage_formal_2,25);
p50wf_2=prctile(wage_formal_2,50);
p75wf_2=prctile(wage_formal_2,75);
p90wf_2=prctile(wage_formal_2,90);


meanwi_2=mean(wage_informal_2);

p10wi_2=prctile(wage_informal_2,10);
p25wi_2=prctile(wage_informal_2,25);
p50wi_2=prctile(wage_informal_2,50);
p75wi_2=prctile(wage_informal_2,75);
p90wi_2=prctile(wage_informal_2,90);



wage_moments=[meanwi_2; meanwf_2;  sdwi_2; sdwf_2;  p10wi_2 ; p10wf_2; ...
    p25wi_2; p25wf_2;  p50wi_2; p50wf_2;  p75wi_2; p75wf_2;...
    p90wi_2; p90wf_2]; % wages

% conditional transitions
formal_2_time1=(state(:,end-1)==2);
informal_2_time1=(state(:,end-1)==3);
nonemp_2_time1=(state(:,end-1)==1);


formal_2_time2=(state(:,end)==2);
informal_2_time2=(state(:,end)==3);
nonemp_2_time2=(state(:,end)==1);


Dfn_2=sum(formal_2_time1==1 & nonemp_2_time2==1)/sum(formal_2_time1==1);
Din_2=sum(informal_2_time1==1 & nonemp_2_time2==1)/sum(informal_2_time1==1);
Dnf_2=sum(nonemp_2_time1==1 & formal_2_time2==1)/sum(nonemp_2_time1==1);
Dni_2=sum(nonemp_2_time1==1 & informal_2_time2==1)/sum(nonemp_2_time1==1);
Dfi_2=sum(formal_2_time1==1 & informal_2_time2==1)/sum(formal_2_time1==1);
Dif_2=sum(informal_2_time1==1 & formal_2_time2==1)/sum(informal_2_time1==1);

transitions=[Dfn_2;Din_2;Dnf_2;Dni_2;Dfi_2;Dif_2];
end





